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The effects of spatial confinements and smooth cutoffs of the waiting time distribution in 
continuous-time random walks (CTRWs) are studied analytically. We also investigate dependences 
of ergodic properties on initial ensembles (i.e., distributions of the first waiting time). Here, we con- 
sider two ensembles: the equilibrium and a typical non-equilibrium ensembles. For both ensembles, 
it is shown that the time-averaged mean square displacement (TAMSD) exhibits a crossover from 
normal to anomalous diffusion due to the spacial confinement and this crossover does not vanish 
even in the long measurement time limit. Moreover, for the non-equilibrium ensemble, we show that 
the probability density function of the diffusion constant of TAMSD follows the transient Mittag- 
Leffler distribution, and that scatter in the TAMSD shows a clear transition from weak ergodicity 
breaking (an irreproducible regime) to ordinary ergodic behavior (a reproducible regime) as the 
measurement time increases. This convergence to ordinary ergodicity requires a long measurement 
time compared to common distributions such as the exponential distribution; in other words, the 
weak ergodicity breaking persists for a long time. In addition, it is shown that, besides the TAMSD, 
a class of observables also exhibits this slow convergence to ergodicity. We also point out that, even 
though the system with the equilibrium initial ensemble shows no aging, its behavior is quite similar 
to that for the non-equilibrium ensemble. 



I. INTRODUCTION 

Recently, slow anomalous diffusion, which is defined by 
the sublinear dependence of the mean square displace- 
ment (MSD) on time, has been found in various phe- 
nomena including lipid granules diffusing in living fission 
yeast cells [1, 2], colloidal particles diffusing on sticky 
surfaces [3] and in networks of entangled actin filaments 
[4], mRNA molecules diffusing in E. coli [5], chromo- 
somal loci diffusing in bacteria [6], telomeres diffusing 
in nuclei of eukaryote cells [7], and proteins diffusing 
in dextran solutions [8]. To understand these phenom- 
ena, two types of slow diffusion models have been exten- 
sively studied so far: (1) continuous-time random walks 
(CTRWs) [1-4] and (2) the generalized Langevin equa- 
tion (GLE) and fractional Brownian motions (FBMs) [5- 
9]. For the CTRW model, if the probability density func- 
tion (PDF) of the waiting times between jumps of the 
particle is a power law w(t) ~ l/r 1+a with < a < 1, 
the ensemble-averaged MSD (EAMSD) shows slow dif- 
fusion [10, 11]: ((Sx) 2 ) ~ t a [12\. Typical origins of 
such power-law waiting times are complex energy land- 
scapes [3, 11, 13-15] and diffusion in inner degrees of 
freedom [16, 17]. Similarly, EAMSD of GLE shows slow 
diffusion if the memory kernel 7/(i) decays algebraically 
rj(t) ~ l/t a with < a < 1 [18, 19]. A possible origin of 
this non-Markovian memory effect is the viscoelasticity 
of the medium [20]. 
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A primary difference between these two stochastic 
processes — CTRWs and GLE — is in their ergodic prop- 
erties. Even if the memory kernel of a GLE is given 
by a power function as stated above, the GLE satisfies 
the ergodic property, that is, time-averaged quantities 
coincide with ensemble-averaged quantities [21, 22]. In 
contrast, this equivalence between the two kinds of aver- 
ages does not hold in the CTRWs with power-law waiting 
times. However, such CTRWs exhibit an extended form 
of ergodicity — time-averaged quantities become random 
variables following a distribution function [15, 23]. This 
distributional ergodicity is called weak ergodicity break- 
ing [15, 23-27] or infinite ergodicity [28, 29] . 

In addition to the abovementioned examples, CTRWs 
are frequently used in various fields of science [10, 11, 30, 
31]. However, there are cases in which some finite-size ef- 
fects should be considered in order to compare the model 
with experimental data. In particular, cutoffs at the tail 
of the power-law waiting time distribution [3, 13, 31] and 
spatial confinement effects [1, 2] arise in many systems. 
The typical origins of waiting time cutoffs are physical 
limits with respect to energy and spatial extensions (e.g., 
[3, 11, 13, 32]). Theoretical analysis of the CTRWs with 
such cutoffs is difficult, because it is necessary to inves- 
tigate the transient behavior. However, this difficulty 
is avoidable by using the tempered stable distribution 
(TSD) P T l(t, A) [22, 23, 33]. The TSD has exponentially 
smooth cutoffs, and is a modified version of the truncated 
stable distribution, which has a sharp cutoff at the tail 
[34]. Furthermore, the TSD has the property called in- 
finite divisibility [22, 33, 35-37], which allows a rigorous 
analysis even for transient behavior. The TSD is a spe- 



2 



cial model for the cutoffs, but it shows typical behavior 
of cutoff distributions such as a slow convergence to the 
Gaussian distribution [22, 33, 35-37]. Furthermore, the 
generalized fractional Fokkcr-Plank equation (GFFPE), 
which we will use in this paper to study EAMSD, has 
been derived for CTRWs with tempered stable waiting 
times [33]. The ergodic properties of the system with- 
out confinement have also been clarified in [23], where 
the existence of a clear transition from weak ergodicity 
breaking (an irreproducible regime) to ordinary ergod- 
icity (a reproducible regime) was shown analytically for 
the case of a non-equilibrium initial ensemble (See Sec. II 
for a precise definition of the non-equilibrium and equi- 
librium ensembles). 

Another important finite-size effect is spatial confine- 
ments. For example, to understand the transport phe- 
nomena in cells [5], confinements due to cell membranes 
should be considered. Confinement effects in CTRWs 
with power-law waiting times were studied numerically in 
[26] and analytically in [38]; it was found that the time- 
averaged MSD (TAMSD) shows a crossover from normal 
diffusion at short timescales to anomalous slow diffusion 
at longer timescales. They also reported numerically the 
weak ergodicity breaking in TAMSD. 

Recently, the CTRW model with the two finite-size ef- 
fects (i.e., waiting time cutoffs and spatial confinements) 
has been used as a model for the transport of lipid gran- 
ules in living fission yeast cells [1, 2], in which confine- 
ment effects are caused by a Hookean force exerted by 
optical tweezers. The model clearly explains the experi- 
mental results such as the crossover in TAMSD and weak 
ergodicity breaking. These studies mainly used numer- 
ical simulations, but a detailed theoretical analysis has 
not yet been reported. Also, they studied only non- 
equilibrium ensemble, and thus the dependences of the 
ergodic properties on initial ensembles are still unknown. 

In this paper, we present theoretical results for CTRWs 
with two finite-size effects: the effects of cutoffs in the 
waiting time distribution and the confinement effects. 
In particular, we focus on TAMSD as an observable, 
and study the crossover from normal to anomalous dif- 
fusion in the TAMSD and ergodic properties in terms 
of the scatter of diffusion constant for the TAMSD. The 
TAMSD, which is often used in single- molecule tracking 
experiments [1, 2, 5, 7, 20, 39-41], is defined as 

J5^ { ^ t) ^_l_j^ A \ x{t ' + A)-x(t'rdt' 1 (1) 

where x{t') is the position of the particle at time t' , t is 
the total measurement time, and A is the time interval. 
Hereinafter we assume that A <C t. He re we define a 
generalized diffusion constant D as (Sx) 2 (A,t) ~ DA@. 
In some experiments, it has been reported that D behaves 
like a random variable depending on each time series [1, 
2, 5, 7, 39]. Therefore, the scatter in TAMSD or D has 
been used to check the consistency of the model with 
experimental data [I, 2, 42]. 



In this paper, we use the TSDs [22, 33, 35-37] as wait- 
ing time distributions of CTRWs. For the TSD, it is 
possible to explicitly write the convoluted waiting dis- 
tributions of any order [Eq. (3)]. Moreover, we use the 
numerical method for the TSD presented in [33] and Ap- 
pendix B, and study the initial and boundary value prob- 
lem of the GFFPE to understand the confinement effect. 

The rest of this paper is organized as follows. In Sec. II, 
we introduce the TSDs. Then, in Sec. Ill, we show the 
crossover from normal to anomalous diffusion in TAMSD 
by using the GFFPE. In Sec. IV, ergodic properties of 
TAMSD are studied using renewal theoretic analysis. 
Sees. V and VI are devoted to conclusion and discus- 
sion. In the appendixes, we summarize some technical 
matters, including a derivation of GFFPE from CTRWs 
as well as the ergodic properties of general observables. 

II. TEMPERED STABLE DISTRIBUTION 

In this paper, we consider CTRWs confined in 
a one-dimensional lattice with unit lattice constant: 
(1,2, ...,L). Jumps are permitted only to the nearest 
neighbor sites without preferences, although this can be 
generalized to jump length distributions with zero mean 
and finite variances. It is also assumed that the par- 
ticle is reflected if it goes beyond the permitted region 
(for example, if the particle jumps into the site L + 1, it 
is pushed back to L). This system is a continuous-time 
version of the discrete-time random walks (DTRWs) with 
reflecting barriers [43]. 

Furthermore, successive waiting times of the particle 
Tfc (k — 1, 2, ...) between jumps are assumed to be mutu- 
ally independent and follow the TSD Ptl {tr , A) defined 
as 

p c\ a —\t 

Ptl(t,A) = -- 

k=i K - 

where Y(x) is the gamma function, a G (0, 1) is the sta- 
ble index, c is a scale factor, and A > is a parame- 
ter characterizing the smooth cutoff. The above defini- 
tion of Ptl (t, A) is a special case of the definition of the 
TSD given in [22, 33, 35-37]. This is because r takes 
only positive values; therefore, we need only one-sided 
distributions. The definition and derivation of TSD is 
presented in Appendix A. When A = 0, Ptl(t, A) is 
the one-sided a-stable distribution with a power-law tail: 
P T l(t,0) - l/r 1+Q as t — )■ oo [44, 45]. Therefore, the 
TSD [Eq. (2)] is the one-sided stable distribution multi- 
plied by the exponential function e _Ar . Thus, Ptl(t, A) 
behaves as Ptl(t, A) <~ e _Ar /r 1+a as r — > oo. Moreover, 
the n-times convoluted PDF P^ l (t, A), which is the PDF 
of the sum of the successive waiting times t n = Y^=i T fc> 
is expressed by using Ptl(t, A) (see Appendix A): 

P£ L (r, A) = n-^PTi/n-^T, n 1/Q A). (3) 
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Thus, the n-times convoluted PDF P^ l (t, A) is also ex- 
plicitly derived from Eqs. (2) and (3). Using Eq. (3), we 
derive transient properties of CTRWs, including various 
crossovers, in the following sections. 

Another advantage of TSD is that its characteristic 
function e^ s ' A ) is also derived explicitly (Appendix A): 

giKCA) =ex p(-c[(A-iC) a -A a ]), (4) 

where e^' A ) = f™ x (1tP tl (t, X)e< T . Equivalently, 

the Laplace transformation of Ptl(t, A), Ptl(s,A) = 
J °° dTe^ ST P T h(T 7 A), is given by 

Ptl(s,A) =cxp(-c[(A + S ) tt -A Q ]). (5) 

Even though we use the TSD Ptl (t, A) as a waiting 
time distribution, we should further specify the first wait- 
ing time Ti, or equivalently, the initial ensemble (Here, 
the first waiting time t\ is the time interval between the 
start of the measurement and the first jump. We always 
assume that measurements start at t' = 0.). In this pa- 
per, we consider two kinds of initial ensembles. The first 
one is a typical and most frequently used non-equilibrium 
ensemble, for which the first waiting times n of the par- 
ticles are chosen from Ptl (t, A) . The second one is the 
equilibrium ensemble, for which the first waiting times T\ 
are chosen from the equilibrium waiting time distribution 
P^ l (t,X). Here, P^(t,X) can be defined by its Laplace 
transformation [46, 47]: 

P^ L ( S ,A)= 1 -^ A) , (6) 

where (r) = c\ a ~ 1 a is the mean waiting time for 
Ptl (t, A) . Note that the second and subsequent wait- 
ing times, T2, T3, . . . , are chosen from Ptl(t, A) for both 
ensembles. Numerical methods to generate random vari- 
ables following Ptl (t, A) and P T £ (r, A) are summarized 
in Appendix B. P^(t,X) can be expressed analytically 
as follows 

e cX " 

x gr^a + i) ( _ cr _ a)fesin(7rfca)/fe(r)) 
fc=l K - 

(7) 

where /fe(r) is defined by /fe(r) = J Q e~ XT l a a ak ~ 1 da. See 
Appendix A for a derivation of Eq. (7). The mean wait- 
ing time for P t \{t,\) is given by (r) = (cA"a + 1 — 
a) I (2A) ~ (1 — a) I (2A), which is much longer than (r) if 
A is small. 

III. ENSEMBLE AVERAGE OF TAMSD 

The confinement effects on CTRWs with power-law 
waiting times were investigated in [38] by using the Frac- 
tional Fokker-Plank equation (FFPE). Here instead of 



the FFPE, we use GFFPE, which was derived in [33], 
to incorporate the smooth cutoff into the waiting time 
distribution and to study the ensemble averages of the 
TAMSD. 

There are three timescales in the present model: (1) 
time interval A, (2) total measurement time t, and (3) 
timescale of the cutoff 1/A. Let us define the Laplace 
variables u and s conjugate to A and t, respectively. 
Since we assume that A <C t (s <C u), it is sufficient 
to consider the following three cases: 

(A) A«t« 1/A A«s«m, 

(B) A « 1/A < t 5«A«n, (8) 

(C) 1/A < A < t s«tt«A. 

Here it is expected that the standard random walk be- 
havior arises in case (C). Thus, we only study the cases 
(A) and (B) in this paper. 

A. Decomposition of ensemble average of TAMSD 

In this su bsecti on, we rew rite th e ensemble averages 
of TAMSD ((Sxf(A,t)) and ({Sxf(A,t)) cq by using the 
EAMSD (\x(A) - x{0)\ 2 ) [see Eqs. (15) and (21) ]. Here 
and in the followings, we use the bracket (•) for the av- 
erage over the non-equilibrium initial ensemble, while 
(■) for the average over the equilibrium initial ensemble. 
Also, we assume that initial position x(0) is uniformly 
distributed on the lattice {1,...,L} for both ensembles. 

First, let w(t) be the waiting time distribution of a 
renewal process with a finite mean, (r) < 00. Moreover, 
we define w e (T;t') as the PDF of the forward recurrence 
time r [47], i.e., w e (r;t') is the waiting time distribu- 
tion at time t' . Here, note that t' is not necessarily a 
renewal time. Particularly, w c (r;0) = w(t) for the non- 
equilibrium initial ensemble, since we assume that we 
start measurements at t' — 0, whereas w e (T;t') 7^ w(t) 
in general. By contrast, for the equilibrium ensemble, 

<V;0 = w cq (r), (9) 

where ui| q (r;t') is the PDF of the forward recurrence 
time at time t' for the equilibrium ensemble, and w cq (r) 
is defined through its Laplace transformation w cq (s) = 
{1 - w(s)}/ (t) s [See Eq. (6)]. This relation [Eq. (9)] is 
obvious because of the time-translation invariance of the 
equilibrium state, and a proof is given in Appendix C. 

1. N on- equilibrium ensemble 

Let us begin with the non-equilibrium ensemble. Using 
the PDF w c (r;i'), we can express the ensemble average 
of the displacement during [t',t' + A] as follows [38]: 

([x(t' + A) - x(t')] 2 ) 

,A 

w / drw c (T;t')(\x(A-r)-x(0)\ 2 ). (10) 
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From Eqs. (1) and (10), we have 

<(fa)*(A,i)>« [ A drw c (T;t)(\x(A-r)-x(0)\ 2 ), (11) 



where w e {r;t) is defined as 



w, 



rt-A 

:(r;i)= / 
Jo 



,«; e (r;i') 



(12) 



From Eq. (11), the ensemble average of TAMSD is given 
by the convolution of w e (r;t) and the EAMSD (|a;(A - 
t) — x(0)\ 2 ). In this and the following subsections, we 
study these two factors using Laplace transformations. 
Here we further rewrite Eqs. (11) and (12) using Laplace 
transformations. First, the Laplace transform of Eq. (12) 
with respect to r gives 



w c (u:t) 



rt-A 

/ dt 
Jo 



, w e (u;t') 
t-A : 



(13) 



where we have defined the Laplace transformations of 
w e (r\t) and w e (r;t) as w e (u;t) = J °° w c (T;t)e~ UT dT 
and w e (u; t) = J °° w c (t; t)e~ UT dr, respectively. Further- 
more, the Laplace transform of Eq. (13) with respect to 
t gives 



w (u; s) = e 



sA 

J s 



ds 



,w e (u; s') 



(14) 



where the double Laplace transformations w c (u;s) and 
w e (u;s) are defined as w c (u;s) = J Q w c (u;t)e~ st dt and 
w c (u]s) = J °° w c (u; t)e~ st dt, respectively. In addition, 
by taking the Laplace transformation of Eq. (11) with 
respect to A, we obtain 



C 



((fe) 2 (A,t))j (u,t) 
= w c (u;t)C[(\x(A)-x(0)\ 2 )] («). 



(15) 



For the non-equilibrium ensemble, the following rela- 
tion between the waiting time distribution w(t) and for- 
ward recurrence time distribution w c (r; t') is well known 
[47, 48]: 



w e (u;s) = 



w(u) — w(s) 



1 



w(s) 



(16) 



See Appendix C for a derivation. If we choose TSD 
[Eq. (2)] for the waiting time distribution w(t), we have 



w e (u; s) ~ 



(a + U ) a - (a + s y 



1 



u — s 



(A + s) a - A" ' 



(17) 



where we have used Eq. (5) and A, it, s -C 1. Using the 
Eq. (17) and s < u, the integral on the RHS of Eq. (14) 
can be approximated as 



as" 



A l-a u a-l 



as 



for A <C s 



for s < A 



(18) 



Then, the inverse Laplace transformation of Eq. (14) with 
respect to s gives 



w e (u;t) ~ < 



r(a + l) 

A l-a u a-l 



t a_1 , for f < 1/A 



(19) 



a 



for 1/A<t 



where we have used e sA ~ 1 [Note that sA <C 1 because 
ofEq. (8)]. 



2. Equilibrium ensemble 

For the case of the equilibrium ensemble, Eq. (10) 
should be replaced by 



([x(t' + A) - x{t')] 2 ) 



eq 



f 

Jo 



dTw cq (T)(\x(A-T)-x(0)\ 2 ). (20) 



Note that the ensemble average in the right-hand side 
(RHS) is taken over the non-equilibrium ensemble. Since 
RHS of the above equation is independent of t', we have 
the following equation through the calculation similar to 
that in the non-equilibrium case: 



C 



((Sx) 2 (A)) cq (u) 



~w^(u)£[(\x(A)-x(0)\ 2 )](u). (21) 

Moreover, the left-hand-side (LHS) of Eq. (20) with t' = 
is just the EAMSD with respect to the equilibrium 
initial ensemble. Thus, we have 

<KA)-.*(0)] 2 ) cq H(^F(A)}c q (22) 

This is a manifestation of ergodicity for the equilibrium 
ensemble. Note that Eq. (22) is valid in general, if (r) is 
finite (and thus w cq (r) exists). 

If we choose TSD [Eqs. (2) and (5)] for the waiting 
time distribution w(t), we have the following relation for 
w c( i(u) [47]: 



w c(1 (u) 



1 



-cl(X+u) a -\ a ] 



(t)u 



a V it 



(23) 



where we used the assumption A -C u [Eq. (8)] and (t) 
cX^a. 



B. Ensemble average of MSD 

Furthermore, we have to calculate the Laplace trans- 
form of EAMSD (|x(A) - x(0)| 2 ) in order to obtain the 
ensemble average of TAMSD [Eqs. (15) and (21)]. First, 
let us express the EAMSD approximately as follows [38] : 

(|a;(A) - a;(0)| 2 ) = t ^ f dxP(x, A; x s , 0)(x-x s ) 2 , 
Jo L Jo 

(24) 
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where P(x, r; x s , 0) is the transition probability from 
x s (t = 0) to x (t = A) of the GFFPE (See Appendix 
E) . It should be note that the GFFPE is an equation for 
the non-equilibrium ensemble, and therefore the EAMSD 
given by Eq. (24) is also for the non-equilibrium ensem- 
ble. The Laplace transformation of Eq. (24) is given by 

C[(\x(A)-x(0)\ 2 )] (u) 



dx s 
o L Jo 



dxP(x, u; x s ,Q)(x — x s ) 2 



6 

I?_ 
6 u 



oo 



7T 



^4 



96 



\u + (nir/L) 2 KuM(u) 

n:odd 



oo 

E 



i 



(A c u) a 



< n 4 (A c u) a + n 2 

n:odd 



(25) 



where we use A< m (A -C 1/A). Moreover, K is defined 
as K = l/2c, and 'n : odd' under the ^2 means that the 
summation is taken over odd terms. We also define a 
characteristic timescale A r as 



A,. 



/2cL 2 



V 7T 2 



1/a 



(26) 



Note that Eq. (25) is the same as the one obtained in 
[38] for the case of the power-law waiting time. This 
means that this property of EAMSD is independent of 
the timescale of the smooth cutoff 1/A. 

Finally, we obtain the following estimate for the 
Laplace transformation of EAMSD: 



for A c u > 1 



C[(\x(A)-x(0)\ 2 )] (u) 



- 2 



— , forA c u<l, 
bu 



(27) 



where we used the RHS of Eq. (25) for the case of A c u <C 
1, while we rewrote the RHS of Eq. (25) by using zeta 
functions (see Appendix F) for the case of A c u 3> 1 as 



1QL 2 



7r 4 A?u 1+Q 



1 



£3 ( A ^) Q 

n:odd 



(28) 




888888 



10 o 
10 



10 



10 10 
Time interval A 



10 



10' 



Q 

I 
< 
H 

o 

<u 
bp 



10' 



io l 



I 10 " 1 



c 



10" 



1 1 

: (b) 


1 1 1 1 




A" 




/o 


o Simulation 






Theory (A«A ) 
Theory (A»A f ) 











10" 



10' 



10 10 
Time interval A 



10' 



10" 



FIG. 1. (Color online) (a) TAMSD (Sx) 2 (A,t) vs. time inter- 
val A in log-log form for the non-equilibrium ensemble (A): 
t < 1/A. The total measurement time t is set as t — 10 6 , and 
the cutoff parameter A as A = 1CP 7 . Other parameters are 
set as a = 0.75, c = 1, and L = 11. TAMSD is calculated 
for 8 different realizations of trajectories and different sym- 
bols correspond to different realizations, (b) The ensemble 
average of TAMSD in log-log form (circles). The lines are 
the theoretical predictions given by Eq. (31). Note that no 
adjustable parameters were used to obtain these theoretical 
lines. 



As shown in the next subsection, the TAMSD for the 
non-equilibrium ensemble behaves differently. Namely, 
the ergodicity is broken even at t — > oo. In contrast, 
the ergodicity is satisfied for the equilibrium ensemble at 
t — > oo. 



C. Ensemble average of TAMSD 



Then, the summation term can be neglected since A c u ^> 
1. From Eq. (27), we obtain the EAMSD for non- 
equilibrium ensemble as follows: 



In this subsection, we derive asymptotic behavior of 
the ensemble-averaged TAMSD using the results from 
the preceding subsections. 



<KA)-*(0)| 2 )^ 



A" 



cT(l + a) ' 
I 6 ' 



for A < A c 



for A > A c . 



(29) 



1. N on- equilibrium ensemble [case (A)]: t <C 1/A 

First, we start with the non-equilibrium ensemble for 
t < 1/A. From Eqs. (15), (19), and (27), we have 
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FIG. 2. (Color online) (a) TAMSD (<5x) 2 (A,t) vs. time inter- 
val A in log-log form for the non-equilibrium ensemble (B): 
t > 1/A. The total measurement time t is set as t = 10 s , and 
the cutoff parameter A as A = 10~ 7 . The other parameters 
are the same as in Fig. 1. TAMSD is calculated for 8 different 
realizations of trajectories and different symbols correspond 
to different realizations, (b) The ensemble average of TAMSD 
in log-log form (circles). The lines are theoretical predictions 
given by Eq. (33). No adjustable parameters were used to 
obtain these theoretical lines. 



slow diffusion at a longer timescale (A c <C A.). As ex- 
pected, these results perfectly coincide with those of the 
previous studies [Eqs. (10) and (11) in [38]]. Also, in this 
regime, the TAMSD depends on the measurement time 
t. That is, the diffusion becomes slower, as the measure- 
ment time increases. We call this behavior aging in this 
article. 

In Fig. 1(a), TAMSDs for 8 different trajectories are 
shown. Although these TAMSDs show similar scal- 
ing b ehavior, the diffusion constant D of each TAMSD 
(Sx) 2 (A,t) ~ DA 7 seems randomly distributed. This 
behavior is analyzed in the next section. In Fig. 1 
(b), the ensemble-averaged TAMSD ((5x) 2 (A,t)) is dis- 
played. The solid and dashed lines are the theoretical 
predictions given by Eq. (31). 



2. N on- equilibrium ensemble [case (B)J: 1/A < < 

Next, we study the case of 1/A <C t. From Eqs. (15), 
(19), and (27), we have 



1 



1 



£ 



((fe) 2 (A,t)) (u,t) = 



1 



L 2 



A" -1 a 6u 2 ~ a 



for A c u > 1 
(32) 

, for A c u < 1. 



Then, taking the inverse Laplace transformation, we ob- 
tain 



cA c 



for A«A C 



((5x) 2 (A,t)) = 



1 



L 2 A l 



(33) 



A Q - 1 a6r(2 - a 



for A c < A. 



£ 



((Sx) 2 (A,t)) (u,t) 

-2 



c r(i + a)t 1 - Q ' 



L 2 



,a-2 



6 r(l + cO^- 



for A c u > 1 



for A c u < 1. 



(30) 



The inverse Laplace transformation with respect to u 
gives 

A 

for A < A c 
(31) 

for A c < A. 



((Sx) 2 (A,t)) 



cT(l + a)t ^ - a, 
L 2 A i- a 

6r(l + a)r(2-a)i 1 - 



Thus, the ensemble-averaged TAMSD shows normal dif- 
fusion at a short timescale (A <C A c ) and anomalous 



Thus, a crossover from the normal to anomalous diffu- 
sion similar to that in case (A) can be observed even 
in the long measurement time limit t — > oo, whereas the 
aging behavior — ^-dependence of the diffusion constant — 
vanishes. Note also that even in this limit, t oo, the 
(ensemble-averaged) TAMSD does not coincide with the 
EAMSD [Eq. (29)]. 

Fig. 2(a) shows TAMSDs for 8 different trajectories. In 
this case, the scatter of TAMSDs, observed in Fig. 1(a), 
is diminished. The ensemble-averaged TAMSD is also 
shown in Fig. 2(b) by circles, where the theoretical pre- 
dictions given by Eq. (33) are shown by solid and dashed 
lines. 



3. Equilibrium ensemble 

For the case of the equilibrium initial ensemble, we 
have the following relation from Eqs. (21), (23), and (27): 
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FIG. 3. (Color online) (a) TAMSD (5x) 2 (A,t) vs. time in- 
terval A in log-log form for the equilibrium ensemble with 
t < 1/A. The total measurement time t is set as t = 10 6 , and 
the cutoff parameter A as A = 10 -7 . The other parameters 
are the same as in Fig. 1. TAMSD is calculated for 8 different 
realizations of trajectories and different symbols correspond 
to different realizations, (b) The same as the figure (a) ex- 
cept that t = 10 8 (t > 1/A). The symbols are TAMSDs for 8 
different trajectories and the thick solid curve is the EAMSD 
(Jx(A) — x(0)] 2 ) . The dotted and dashed lines are theoret- 
ical predictions given by Eq. (35). Note that no adjustable 
parameters were used to obtain these theoretical lines. 



The equation (35) is equivalent to Eq. (33), the 
TAMSD for the non-equilibrium ensemble [case (B)]. 
However, Eq. (35) is valid for arbitrary measurement 
times t, while Eq. (33) is valid only for long measure- 
ment times (1/A -C t). In addition, the aging behavior is 
absent in the equilibrium case. 

In Fig. 3(a) and (b), TAMSDs for 8 different trajecto- 
ries are shown for a short and long measurement times t, 
respectively. Surprisingly, the scatter of TAMSD is even 
broader than that in the non-equilibrium case [Fig. 1(a)] 
at short measurement times t as shown in Fig. 3(a). 

In summary, for the non-equilibrium ensemble, the 
scatter of TAMSD appears with the aging behavior 
[the non-equilibrium ensemble (case A)] , whereas for the 
equilibrium ensemble the scatter appears without aging. 
These scatters are also studied quantitatively through 
theoretical analysis in the next section. The theoretical 
predictions given by Eq. (35) are also shown by dotted 
and dashed lines in Fig. 3(b). 

Furthermore, from Eq. (22), the EAMSD for the equi- 
librium initial ensemble (jx(A) — x(0)] 2 ) c is also given 

by the RHS of Eq. (35). In Fig. 3(b), a numerically 
obtained EAMSD ([x(A) - x(0)] 2 ) cq is displayed by a 

solid curve, which is consistent with the theory [Eq. (35): 
the dotted and dashed lines in Fig. 3(b)]. Thus, the 
EAMSD for the equilibrium ensemble ([x(A) — x(0)] 2 ) eq 

[Eq. (35)] is different from the EAMSD for the non- 
equilibrium ensemble ([x(A) — x(0)] 2 ) [Eq. (29)]. 



IV. STATISTICAL AND ERGODIC 
PROPERTIES OF TAMSD 



The statistical property of TAMSD is dominated by 
the property of the number of jumps Nt until time t (see 
Sec. IV C). Therefore, we first study N t in Sees. IV A 
and IV B on the basis of the analysis presented in [23] . 



for A c u > 1 



C 



((Sx) 2 (A)) cq (u) 



1 1 

c\ a ~ 1 a u 2 ' 

iZ=T~ a 2-a ' fOT AcU <C L 



(34) 



Through the inverse transformation, we have 



A 



<(&0 2 (A)>« 



1 L 2 A X -" 
. \ a - 1 a 6r(2 - a) 



for A < A c 
, for A > A c 



(35) 



A. Real space analysis 

In the Sec. II, we define t n as the time when the n-th 
jump occurs for a trajectory x(t): t n = X)fc=i T fc- From 
this definition, we have the following relation: 

G(n; t) = Prob (N t < n) 

= Pror(t„ > t) = Prob ^ T k > tj , (36) 

where Prob(...) is the probability and Tf. (k — 1, 2, ...) are 
the waiting times between jumps. 
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1. N on- equilibrium ensemble 

From Eq. (36), we have the following equation for the 
non-equilibrium ensemble: 

/oo 
O 

dr F T L(r,n 1/Q A), (37) 

-i/» t 

where we used the statistical independence between the 
waiting times {k = 1,2, ...) and Eq. (3). Furthermore, 
if we change variables from n to x as n = t a x, we obtain 

Prob <x) = J dr P TL (r,tx^ a \j 

(38) 

Note that the integrand of the RHS of Eq. (38) is not the 
PDF, because it contains the variable x. To derive the 
PDF for x, we insert Eq. (2) into Eq. (38), and then we 
have 



Pm!> ( — \ < x 1 = 



rv7r — ' 



fc=l 



x (— cx) k sin(7r ka)ak, 



where 



ak 



dre 



(39) 



(40) 



Differentiating Eq. (39) in terms of x, we obtain the PDF 
for x: 



e^f r(fca + l). 

h(x; t) = — - fc , H -c)" 

fe=i 



CV7T 



c(U) a a; 



aT 1 sin(7rfco;)afe. (41) 



When A = 0, this PDF f (x) = f {x,t) is called the 
Mittag-Leffler distribution [28, 29]. 

As shown below [Sec. IV C], TAMSD behaves as 
(fe) 2 (A, t) « ^A for small A. Thus, the diffusion con- 
stant D t , defined as {Sx) 2 {A,t) DtA, shows the same 
statistical property as that of N t . Fig. 4 shows the PDF 
of D t for two different measurement times t. Note that 
the PDF narrows as t increases. The analytical results 
given in Eq. (41) are also depicted by solid and dashed 
lines in the figure. 



2. Equilibrium ensemble 

The equilibrium case is also analyzed in a similar way. 
In this case, Eq. (37) should be replaced by 



Prob {N t < n) 



dr {P^l*P^ X ){r,X), (42) 




1 1.5 
Diffusion coefficient D t 




1 1.5 
Diffusion coefficient D, 

FIG. 4. (color online) (a) PDF of the diffusion coefficient Dt 
for the non-equilibrium TAMSD {5x) 2 {A, t) fa D t A for small 
A. Each PDF is normalized so that its mean value equals 
unity. D t is calculated from TAMSD by least-square fitting 
over the interval < A < 10. The results for two differ- 
ent values of measurement times are presented: t = 3 x 10 6 
(circles) and 3 x 10 7 (squares). The other parameters are 
set as A = W~ 7 ,a = 0.75, c = 1, and L — 11. The lines 
correspond to the theoretical predictions given by Eq. (41). 
No adjustable parameters were used to obtain these curves, 
(b) PDF of the diffusion coefficient D t for the equilibrium 
TAMSD {Sx) 2 {A,t) w A A for small A. The parameter val- 
ues are the same as in the figure (a). The lines are to guide 
the eye. 



where {f*g) means a convolution. In contrast to the non- 
equilibrium case, however, it seems difficult to obtain a 
simple expression of the PDF of Nt for the equilibrium 
case. As shown in Fig. 4, a qualitative difference ap- 
pears for short measurement time regime, t < 1/A. In 
fact, there is a peak at D t — for the equilibrium ini- 
tial ensemble, which is due to the mean waiting time of 
P£l(t, A), (t) (~ 1/A), is much longer than that of 
Ptl(t, A), (r) (~ 1/A 1_Q ). Namely, there are many tra- 
jectories that are trapped more than the measurement 
time t {< 1/A) [49]. 
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B. Laplace space analysis 

Next, to clarify the ergodic properties of the sys- 
tem, we study the relative standard deviation (RSD) 
R(t) = y/(N?)J (N t ), where (-) c is the cumulant. The 
quantity R(t) is a measure of ergodicity [26]. In fact, if 
R(t) — 0, time averages of an observable give the same 
value independent of the trajectory; however, if R(t) > 0, 
they are different from one trajectory to another. This 
quantity R(t) was also used in some molecular dynamics 
simulations to characterize the ergodicity breaking and 
Gaussian fluctuations of lipid motions in cell membranes 
[50] , and to determine the longest relaxation time in en- 
tangled polymers [50]. 



1. N on- equilibrium ensemble 

Let us begin with the non-equilibrium ensemble. In 
order to derive an analytical expression of R{t) , we start 
with the Laplace transformation of Eq. (37): 



G(n;s) = / dte- ts 
Jo 



f-OO 

Jn-i/° 



dr P Th (T,n^ a X) 



1 _ e -„ c [(A+s) a -A Q 



(43) 



where G(n; s) is the Laplace transformation of G(n; t) 
in terms of t. Next, let us define a function g(n; s) as 
g(n; s) = G(n + 1; s) — G(n; s). Note that g(n; s) is the 
Laplace transformation (with respect to t) of the PDF of 
N t . Furthermore, we perform a discrete Laplace trans- 
formation of g(n; s) with respect to n as follows: 

- q[v . :) - 1 l-axp(-c[(A + »)«-*■]) {u) 
9[ ' ' ~ s 1 - exp {-v - c[(A + *y - A"]) ' [ ' 

where we define g(v\ s) as g(y\s) = X)^Lo e ~ nv g{ n \ s )- 
Using the assumption s, A, v <C 1, we have 



1 00 k 



(45) 



fc=0 



a. First moments From Eq. (45), we obtain the 
Laplace transform C[{N t )]{s) of the first moment (Nt) 
as 



£[(Nt)](s) 



1 



cs 



a+1 ' 
1 



s > A 



c\ a ~ 1 as 2 



1 + (1 - «)^J , s< A. 



(46) 



The inverse Laplace transform of the above equation is 
given by 



t° 



(Nt) 



cT(a + iy 

t 1 — a 



cX^a 2cX a a 



t <1/A 
t >1/A. 



(47) 



The ensemble-averaged mean square displacement 
(EAMSD) for free CTRWs (i.e., CTRWs without con- 
finement effects) is known to be proportional to (Nt) [11]: 
(((5r) 2 ) (t) - (N t ). Thus, the EAMSD of the present 
model without the effect of confinement shows transient 
subdiffusion, i.e., subdiffusion for short timescales and 
normal diffusion for long timescales [13]. The crossover 
time between these two regimes is characterized by 1/A. 

b. Second moments Similarly, the Laplace trans- 
form £[(iV t 2 )](s) of the second moment (-/V t 2 ) is given 

by 



c 2 s 2a+l ' 

2 

V. C 2^2a-2Q,2 S 3 



l + (l-a)- 



s > A 

, s < A. 
(48) 



Then, the inverse transform is given by 
2t 2a 



c 2 r(2a + 1) 
2 

I c 2 A 2a " 2 a 2 



t < 1/A 



t 2 I- a 



(49) 



t > 1/A. 



c. Relative standard deviation Using Eqs. (47) and 
(49), we obtain asymptotics of the RSD y/(N?)J (Nt): 



V<MT C 

(Nt) 



~ < 



'2r 2 (a + l) 



r(2a + l) 



'l-a 
Xt 



1, t<l/A 



t > 1/A. 



(50) 



We define a crossover time t c between the two regimes, 
t <C 1/A and t >• 1/A, as the intersection of the two 
functions in Eqs. (50); therefore, we have 



t c = 



(1 



2r 2 (q + l) 



r(2a+l) 



1 



-A" 



(51) 



As shown in Fig. 5, the RSD remains almost constant 
before the crossover time t c , and starts decaying rapidly 
after the crossover. In the figure, the RSD for the ex- 
ponential waiting-time distribution, which has the same 
mean waiting time (t) as that of the TSD with A = 10~ 6 , 
is also shown by pluses. It is clear that the RSD for 
the exponential distribution (pluses) decays much more 
rapidly than that for the TSD (triangles). 



2. Equilibrium ensemble 

The calculation for the equilibrium case is almost par- 
allel to that for the non-equilibrium case except that we 
should use Eq. (42) instead of Eq. (37). Thus, we only 
show the final results. First, the generating function 
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FIG. 5. (color online) (a) RSD y/{D?)J (D t ) vs. total mea- 
surement time t for the non-eq uilibriu m initial ensemble. Dt 
is calculated from the TAMSD (Sx) 2 (A,t) by a least-square 
fitting over the interval < A < 1. Three different values 
of A are used: A = 10 -4 (circles), 10 -5 (squares), and 10 -6 
(triangles). In addition, a, c, and L are set as a = 0.75, 
c — 1, and L — 11, respectively. The lines correspond to the 
theoretical predictions given by Eq. (50); the solid line is the 
result for short timescales t < i c , and the dashed lines are 
for long timescales t S> t c . The intersections of the solid and 
dashed lines correspond to the crossover times t c given by 

Eq. (51). (b) RSD ^ {D% ) eq J (A) eq vs. total measurement 

time t for the equilibrium initial ensemble. The parameter 
values are the same as those in the figure (a), except that the 
results only for two different values of A are plotted for clar- 
ity: A = 10 -5 (squares), and 10 -6 (triangles). The lines corre- 
spond to the theoretical predictions given by Eq. (54). (a) and 
(b) The plus signs in both figures are the RSD for the case in 
which the waiting time distribution is given by the exponen- 
tial distribution P(t) — exp(r/ (t))/ (t) with the same mean 
waiting time as that of the TSD with A = 10 -6 (triangles): 
(r) = c\ a ~ 1 a. The dot-dashed line is a theoretical predic- 
tion for the exponential distribution: R(t) = (caA a ~ 1 /i) 1/ ' 2 . 
Note that the scales of vertical axes of two figures are differ- 
ent. For comparison, the theoretical prediction of RSD for 
the non-equilibrium ensemble at short timescales is depicted 
by a dotted line in the figure (b). 



{v; s) is given by 



i-fc+1 



Using this function, we obtain (N(t)) 



cq 



t 



cq cA° 



(52) 
(53) 



and the RSD for the equilibrium initial ensemble 



W 2 >e q ,c 



2a 



T(2 + a)(Ai) 



(Nt. 



(54) 



'1-q 
Xt 



t> 1/A. 



Although the RSD slowly decays at the short measure- 
ment timescales, the relative fluctuations are even larger 
than the non-equilibrium case as shown in Fig. 5. 



C. Statistical properties of TAMSD 

In this section, we discuss the statistical properties of 
TAMSD for one-dimensional CTRWs, although the anal- 
ysis for free CTRWs in the followings can be generalized 
for higher dimensional systems [23]. Properties applica- 
ble to a general class of observables are summarized in 
Appendix G. 

To begin with, the TAMSD can be approximately given 
by the time average of the following observable: 



h(t') = J25(t' -t k )h k 



(55) 



k=l 



fe-1 



h k = AzZ + 2^2z k zt6(A-(t k -tt)), (56) 



i=i 



where z k — ±1 is the displacement of the jump at time 
t' = t k , and 9(t) is a step function defined by 



0(t) 



0, (t < 0) 
t, (t > 0). 



(57) 



Here we prove that the time average of hit') defined by 
Eq. (55) is the TAMSD. Let us express a trajectory of a 
CTRW as 



;(*') = J2 ZkI ( tk < 0. 



(58) 



k=l 



where I(t k < t') is defined as I(t k < t') = 1 if the inside 
of the bracket is satisfied, while I(t k < t') = otherwise. 
Then, the displacement x(t' + A) — x{t') is given by 



:(*' + A) - x(t') = ZkHt' <t k <t' + A). (59) 



fc=i 
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Furthermore, the squared displacement [x(t'+A)— x(t')} 2 
is expressed as 

oo 

[x(t' + A) - x{t')] 2 = Kh -A<t'< t k )z 2 k 



k=\ 



oo k— 1 



+2j2J2 ZkZiI( - tk - A<t ' <ti ^ 

(60) 



fc=l 1=1 



From Eqs. (60) and (1), wc obtain the following approx- 
imation for the TAMSD: 



(5x)Z(A,t)~-J2 



k=l 



k-1 



Azl + 2^2z k zi0(A-{t k -t t )) 



(61) 

Now, it is clear that the RHS of Eq. (61) is equivalent to 
the time average of h(t'), defined by Eqs. (55) and (56). 



1. CTRWs without reflecting boundaries 

In the absence of the confinement effect [23] (i.e., with- 
out reflecting boundaries), the relations (hk) = A and 
(hkhk+n) — (hk) (hk+n) = for n > 1 hold because of the 
mutual independence of z k , (z k ) = and z\ = 1. In ad- 
dition, since z k is independent of the initial ensemble, the 
same relations can be obtain also for the equilibrium en- 
semble: (Mcq = A aild ( h kh k +n) cq - (hk) cq (hk+n) cq = 

0. Then, h k [Eq. (56)] satisfies the law of large numbers 
[Eq .(G2)]; therefore, we have 



i N * n ^ Nt n 

(fe)2(A,t) Rj -^/ lfe = ^-^/ lfe = ^A, (62) 



k=i 



fe=i 



for arbitrary A. Note that this result is independent of 
the choice of the initial ensembles. Therefore, TAMSD 
increases in proportion to A, and has the same statistical 
property as N t . In particular, the diffusion constant D t 
is given by 



D t = 



(63) 



Furthermore, the PDF of D t is given by Eq. (41), and 
the RSD of D t shows the asymptotics such as the ones 
given in Eq. (50). 



2. CTRWs with reflecting boundaries 

Although a similar property holds for the system with 
reflecting boundaries, the calculation becomes more com- 
plicated, because z k {k = 1, 2, ...) are not mutually inde- 
pendent. Here, wc present only an outline of the proof. 
The main tool is the spectral decomposition of the n- 
time transition matrix P^> of the DTRWs with reflecting 



boundaries [43]: 



L-l 



f*0 = £| r ) {r \ K . 



(64) 



r=0 



where L is the total number of sites. The ket |r) is the L- 
dimensional eigenvector of the transition matrix defined 

by 

/ L l/2 \ 



|o> 



(65) 



for r = 0, and 
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ri 



\ r L, 



with 



sin ■ 



sin ■ 



LV2(1 



cos ■ 



(66) 



for 1 < r < L — 1. Here, j is a site index (1 < j < L). 
Also, A r are the associated eigenvalues: Xq = 1 for 
r = and A r = cos(nr/L) for 1 < r < L — 1. Thus, 
the first term (r = 0) in Eq. (64) corresponds to the 
eigenmode with the unit eigenvalue (the non-decaying 
mode), which is a uniform state, and the other terms 
decay exponentially fast to zero. Each element of the 

(n) 

transition matrix, i.e., • , is the transition probability 
from the site i to the site j during n jumps. An impor- 
tant point is that each element of the matrix tends to 
1/L exponentially fast as n increases. It follows that 
correlation functions for the DTRWs decay exponen- 
tially, too. By using this fact, we can show (after some 
lengthy calculations) that the correlation functions for 
the CTRWs (hkhk+ n ) — (hk) (hk+n), which can be ex- 
pressed with correlation functions for the DTRWs such 
as (z n ziZkZ\) — (z n zi) (zkZ\), also decay exponentially. 

As a result, the law of large numbers holds for the hk 
even for the confined system, and we obtain 

AT 1 N * AT 

(fe) 2 (A, t) « -f— hk = -fun (A) . (67) 

* k=l 

The function /Ufc(A) can be derived by using the tran- 
sition matrix [Eq. (64)] with a suitable hydrodynamic 
limit, however an easier way here is to utilize the results 
in Sec. III. For example, from Eqs. (31) and (47) [or 
Eqs. (33) and (47)] we obtain 



($x) 2 (A,t) 



N t cI?A x - a 
lT6r(2-a)' 



for A < A c 



for A > A c 



(68) 



for the non-equilibrium initial ensemble. Similarly, from 
Eqs. (35) and (53), we have Eq. (68) also for the equi- 
librium ensemble. Thus, Eq. (68) is valid for both non- 
equilibrium and equilibrium ensembles. Thus, all the dif- 
ferences between the two ensembles are included in the 
statistical properties of N t in Eq. (68). See, for example, 
Eqs. (47) and (53). 
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V. CONCLUSION 

Up to now, two characteristic properties have been re- 
ported for confined CTRWs [38]: (i) there is a crossover 
from normal to anomalous diffusion, and (ii) TAMSDs 
are distributed depending on trajectories (i.e., weak 
ergodicity breaking). These results are for the non- 
equilibrium ensemble [There is no equilibrium ensemble 
for CTRWs with power-law waiting times with stable in- 
dex a £ (0, 1), because the system never reaches an equi- 
librium state.]. In this paper, in addition to this con- 
finement effect, we incorporated a cutoff into the waiting 
time distribution. For the case of the non-equilibrium 
ensemble, we analytically determined that property (i) 
persists even in the long measurement time limit t — > oo 
(This fact was first found numerically in [2]). In con- 
trast, the property (ii), that is, the distributional behav- 
ior of TAMSDs, appears for short measurement times t, 
whereas ergodicity holds for longer measurement times. 
The important point is that as compared to common dis- 
tributions such as exponential distribution, it takes a very 
long time to observe ergodic behavior for the case where 
the waiting time distribution is given by the TSD. In ad- 
dition, this transition from weak ergodicity breaking to 
ergodicity is a transition from an irreproducible regime 
to a reproducible regime. 

Furthermore, as shown in Eq. (51), the crossover 
time t c between weak ergodicity breaking and ergodic- 
ity is proportional to 1/A [also, the crossover time be- 
tween anomalous and normal diffusion in EAMSD for 
free CTRWs is proportional to 1/A; see the text below 
Eq. (47)]. Because the mean waiting time (r) is given by 
(t) - I/A 1 "", we have 

*c~<t>^. (69) 

Thus, the crossover time t c is not proportional to the 
mean waiting time (t) ; in fact, t c can be much longer 
than (r) for a close to 1. These facts may be important 
for estimating the crossover time in experiments [3]. 

In contrast to the CTRWs with power-law waiting 
times with stable index a € (0, 1), the CTRWs with cut- 
off waiting times have an equilibrium state. It is surpris- 
ing that the crossover from normal to anomalous diffusion 
[property (i)], and a scatter in TAMSD [property (ii)] ex- 
ist even for the case of the equilibrium initial ensemble. 
The scatter in TAMSD is even broader than the case 
of the non-equilibrium ensemble. The main difference 
from the non-equilibrium case is that there is no aging 
for the equilibrium ensemble. Another difference between 
the two ensembles is that the RSD decays algebraically 
even in the short time regime for the equilibrium case, 
whereas it does not decay for the non-equilibrium case. 
These properties for the equilibrium ensemble might fit 
well with some experimental results [1, 2, 5, 7]. The 
important point is that the absence of aging in experi- 
mental data does not necessarily exclude the possibility 
of the CTRWs. We also presented a numerical method 
to generate the equilibrium ensemble (Appendix B). 



VI. DISCUSSION 

As shown in Appendix G, time averages for a class of 
observablcs including TAMSD follow the ML distribu- 
tion. However, different distributions can arise for differ- 
ent classes of observables [51]. As an example, we show in 
Appendix G the case in which long time averages follow 
the PDF called the bilateral ML distribution [52]. 

Recently, anomalous properties have been found in 
both experiments [1-3, 5-8, 53] and molecular dynam- 
ics simulations [50, 54]. The CTRWs, which are a model 
of random walks in random environments, are often used 
as a model of these anomalous properties. However, there 
are qualitatively different types of random walk models 
in random environments: for example, barrier models 
[11, 55], random force models [11, 56] and the reptation 
model for entangled polymers [50, 57]. To the best of the 
authors' knowledge, the ergodic properties, such as the 
behavior of RSD, of these systems are still unclear. 

Besides random walks in random environments, there 
are still several different mechanisms for anomalous be- 
havior, e.g., GLE, FBMs, and diffusion on fractal struc- 
tures [58]. Therefore, it is important to develop tech- 
niques for time series analysis to elucidate which mech- 
anism is the actual cause of the anomalous properties 
observed in various experiments and molecular dynami- 
cal simulations [9, 59]. It is also important to investigate 
systems in which these mechanisms are combined [60], 
and those in which nonlinear dynamics are coupled with 
these anomalous mechanisms [21, 53]. 

Furthermore, the TSD used in this article has the infi- 
nite divisibility, which seems to be related to a renormal- 
ization group transformation [36]. Thus, it is expected 
that the TSD has a kind of stability (or instability) in a 
functional space of probability densities. It is important 
to clarify such stability of the TSD. 
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Appendix A: Infinitely divisible distribution 

Here we briefly introduce the one-sided TSD. In gen- 
eral, the infinitely divisible distributions P(t) is defined 
as follows: if oj(s) is the characteristic function (Fourier 
transform) of P(t), then there exists a characteristic 
function cj„(s) such that to(s) = to™(s). Note that this is 
not trivial, because the function uj 1 ' n {s) is not necessar- 
ily a characteristic function of a PDF (i.e., a non-negative 
function) . 



13 



Now let us define the one-sided TSD Ptl(t, A). First, 
we define the characteristic function e^^' X > of Ptl{ t , A) 
as follows: 



TL 



1 f c 



V(C,A)„-iCr 



e-^dC, (Al) 



oV>(C,A) 



/oo 
P T L(T,A)e lCT dT. 
-oo 



(A2) 



Then, the function i[)(s, A) is defined by the canonical 
form of the infinitely divisible distributions [61]: 

/oo 
(e^-l)f(r,X)dr. (A3) 
-oo 

If /(r, A) is a PDF in terms of r, then e^> A ) is the char- 
acteristic function of the compound Poisson distribution 
[61]. For the TSD, however, the function /(r, A) is not a 
PDF, and is defined as follows [22, 33, 35-37]: 



f(r,X) = { r- 1 "^' 



T(-a) 



(r<0) 
(r > 0), 



(A4) 



where c is the scale factor, and a and A are parame- 
ters satisfying < a < 1 and A > 0, respectively. By 
changing the integration path on the complex plane in 
Eq. (A3) to C : z = s(X + i()/(X 2 + C 2 ) [fl € [0, oo)], we 
obtain A) as 

lKC,A) = -c[(A-*C) a -A Q ]. (A5) 

The above equation indicates an important property: 

nV(C,A)=V(n 1 /%n 1 / a A), (A6) 

which is related to the infinite divisibility introduced 
in the first paragraph of this section, because it can 
be rewritten as ij)( n ~ 1 ^ a Ci n~ x / a \) — ip((,X)/n, and 
?/>(n~ 1//Q C, n~ 1 / a X) is obviously a characteristic function. 
It is clear that Eq. (3) is obtained from Eq. (A6). 

Now, we derive Eq. (2). From Eqs. (Al) and (A5), we 
have 



Ptl(t,A) = — Re j « 



„cA° 



Im / e-< x -^ a e- s ds, (A7) 
Jo 



where we have changed the integration path from the 
real to the imaginary axis by using the Cauchy's integral 
theorem. Then, we have 



Ptl(t,A) 



7TT 

cA a — At 

7TT 



Im / 

tA 

oo 

Im / e 



e s ds 

rA 

/ e- c (--)"e- s ds. (A8) 
Jo 

Finally, we obtain Eq. (2) by Taylor expansion of the 

exponential function e and the integral represen- 

tation of the Gamma function T(ka + 1) = / °° s afc e~ s . 



The equilibrium waiting time distribution (7) can be 
derived by using Eq. (B5) as follows: 



Jo o,{t) 



iL[ -,X)da. (A9) 



Appendix B: Method of Numerical Simulation 

Let Y\ be a random variable following the TSD 
[Eq. (2)]. We briefly review the method to generate Y\ 
numerically according to [33]. First, a random variable 
Yo, which follows the one-sided stable distributions, can 
be generated by the following equation [62]: 



Y 



Ma 



sin (a 

(cos!0 1/Q 

cos(V-a(V+f)) 
W 



(l-a)/a 



(Bl) 



where V is a uniform noise in [— 7r/2, 7r/2], and W is an 
exponential noise with mean 1. To generate Y\, note that 
TSD [Eq. (2)] is just a stable distribution multiplied by 
the exponential factor e~ Xr . Therefore, we first generate 
Yq, then we accept it with a probability exp(— XY ) < 
1. If Yo is rejected, it will be regenerated according to 
Eq. (Bl) until it is accepted. Thereafter, the random 
variable that is finally accepted follows TSD [33]. 

Similarly, to generate equilibrium noise Y^ 9 , first we 
create an auxiliary random variable T by accepting Yq 
with a probability AYb exp(— XY + 1). Then, Y^ 9 is ob- 
tained by Y^ q = XT, where X is a uniform noise in [0, 1]. 
This can be proved as follows. 

First, let us consider a time axis (—00,00) covered by 
many non-overlapping time intervals. Each interval is ob- 
tained from Ptl (t, A) (The intervals are assumed to be 
mutually independent). Then, choose a time t randomly 
from the time axis. Now, the random variable T is de- 
fined by the length of the time interval which includes 
the randomly chosen time t. In addition, let p(r) be the 
PDF of the variable T. Then, p(r) is given by 

p(r) = ^-Ptl(t, A) cx re- XT P Th (T, 0), (B2) 

because the probability with which an interval is cho- 
sen is proportional to the length of that interval r. 
Thus, the random variable T can be generated from 
the stable noise Yo by accepting it with the probabil- 
ity XYo cxp(— AY" + 1) < 1 (The factor Ae is added to 
make the maximum probability be unity for computa- 
tional efficiency). Moreover, by the Laplace transform 
of the above equation, we have 



(B3) 
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Next, let us define a new variable Y^ q = XT, where X 
is the uniform noise in [0, 1]. Then, we have 

Prob (F A oq > r) = Prob (XT > r) = f Prob (t > -) da. 

Jo ^ a ' 

(B4) 

By differentiating the both side, we obtain 



Pcq(r) 



r\ da 



(B5) 



where p e q(T~) is the PDF of Y^ q . The Laplace transform 
of the above equation leads to 

Pcq(s) = / p(as)da. (B6) 
Jo 

Inserting Eq. (B3) into the above equation, we obtain 
Peq(s) = [1 — Ptl(s,X)]/ (t) s, which is the Laplace 
transform of the equilibrium waiting time PDF [Eq. (6)] . 
Therefore, Y^ q follows the equilibrium waiting time dis- 
tribution P,£l (t, A) . 

This method to generate the equilibrium waiting time 
might be applicable to other probability densities with 
finite moments. 



Appendix C: Time translation invariance of the 
equilibrium waiting time 

Let w(t) be a waiting time distribution for a renewal 
process, and w eq (r) be the associated equilibrium waiting 
time distribution. Here, we show the relation u>| q (r; t') — 
w cq (r) [Eq. (9)] by the method presented in [47]. First, 
let us define 

w(t; t',n) = (S (r - (t n+1 - t')) I [t n < t' < t n+1 )) , 

(CI) 

where I(a < t' < b) = 1 if the inside of the bracket is 
satisfied, otherwise. In addition, t n (n — 1,2,...) arc 
renewal times defined by t n = Y^k=i T k i where are 
successive waiting times. The Laplace transform with 
respect to r and t' gives 



uitr.x.n) = ( J e -( t »+ 1 - t '>"e- t ' a dt / 




for n = 

i 

^(^(.J.MlM, forn>l. 
s — u 

(C2) 



Because w(u; s) — X^^Lo s ' n )' we nave 

x w(u) — w(s) w cq (s) w cq (u) — w eq (s) 
ffl">; s) = — — — • rr^r + 



s — u 1 — w(s) 



s — u 



Here, if we replace w eq (s) with w(s), we obtain the for- 
ward recurrence time distribution for the non-equilibrium 
ensemble [Eq. (16)]: 



w c (u; s) 



w(u) — w(s) 1 



s — ii 1 — w(s) ' 



(C4) 



It follows that 



w eq (u) = lim sw c (u;s) = . ( C 5) 

s^O (t) u 

By inserting Eq. (C5) into Eq. (C3), we find that 
u>l q (u;s) — w eg (u)/s. The inverse Laplace transform 
gives the relation we want. 



Appendix D: Derivation of GFFPE 

GFFPE was derived from a subordinated process in 
[33]. Here, we derive GFFPE from CTRWs in a hydro- 
dynamic limit [11, 63, 64]. Let P(x,t) be the PDF of 
a particle at a time t, and Po(x) be the initial density 
Po(x) = P(x,0). (In this section, we consider t to be 
the usual time variable, instead of the total measurement 
time). Moreover, each particle is assumed to follow the 
CTRW dynamics on the real line x or a one-dimensional 
lattice x — xq,x±i,x±2, .... Now, we define ip(x,t)dxdt 
as the probability of a particle to perform a jump with 
length x after being trapped for a certain period t. Then, 
the probability density of a particle to be trapped for a 
period t is given by 



dx' ip(x',t')dt' = 1 - / w{t')dt'.(m) 

-oo JO Jo 

Furthermore, we define Q(x,t)dtdx as the probability of 
a particle to reach an interval [x,x + dx) in the period 
[t, t + dt). Then, we have 

P(x,t)= f dt'(j>{t-t')Q(x,t') + cj>(t)P (x), (D2) 
Jo 

/oo />£ 
dx' / dt'ip(x' ,t')Q(x- x' ,t-t') 
-oo Jo 

dx'i)(x',t)P {x-x') (D3) 



Taking the Fourier and Laplace transforms with respect 
to space and time, respectively, we obtain 



Q(k,u) = 



P(k,u) 



4>(k,u)P (k) 

1 — tp(k, u) 
l-w(u) P (fc) 



(D4) 



(D5) 



(C3) 



where we used the relation 4>(u) = (1 — w(u))/u [the 
Laplace transform of Eq. (Dl)]. 



15 



Next, we assume that the PDF tp(x,t) can be separa- 
ble, i.e., tp{x,t) = l(x)w(t). We further assume that 



l(k) ~ 1 - 



(D6) 



in the hydrodynamic limit k — > 0, where (fa 2 ) is 
the mean displacement of a single jump, i.e., (fa 2 ) = 
dx x 2 l(x). For one-dimensional CTRWs with jumps 
only to the nearest neighbor sites, we have (fa 2 ) = 1. 
Moreover, we use the TSD, -Ptl(£) ; as the waiting time 
distribution w(t): 

w(u) ~l-c[(\ + u) a - A"], (D7) 

where A, u <C 1 is assumed. Under these assumptions, 
Eq. (D5) can be rewritten as follows: 



uP{k,u) - P (k) ~ - 



(fa 2 ) 



uk 2 



2c [(\ + u) a -\ c 



rP(fc,ti).(D8) 



The inverse Fourier and Laplace transformations with 
respect to space and time lead to 



dP(x,t) =K ^dP(x,t) 



dt 



d 2 a 



(D9) 



where K is a constant given by K = (fa) /2c, and <l t is 
an operator defined by 

$ t f(t) = j t J*M(t-t')f(t')dt' (D10) 

with function M(t), which is defined by its Laplace trans- 
form: 



M(u) 



dte- ut M(t)dt = 1 ;-.(Dll) 

w (A + u) a - X a ' 



As A -> 0, Eq. (D9) leads to the usual FFPE [38, 64]. 
Thus, the above equations (D9)- (Dll) is called gener- 
alized FFPE [33]. Note also that GFFPE is a equation 
for the non-equilibrium initial ensemble. In contrast, we 
do not know what kind of equation describes the system 
started with the equilibrium ensemble. 



Appendix E: Mixed problem of GFFPE 

The initial-boundary value problem of GFFPE 
[Eq. (D9)] under the boundary condition 

dP , N dP , N 

- M = -(L,t)=0, (El) 

and initial condition 

P(x,0) = 6(x-x s ) (E2) 

can be solved by a standard method for the diffusion 
equation [38]. In fact, by assuming the separability 
P(x,t) = X(x)T(i), we obtain 



T'{t) _ X"(x) 



KQ t T(t) X{x) 



= -q, 



(E3) 



where q is a constant. X(x) satisfying the boundary con- 
dition [Eq. (El)] is given by X n (x) = cos(mrx/L) and 
q = (mr/L) 2 (n = 0,1,2...). Similarly, T(t) is derived in 

T n (0) 



its Laplace form: 

f n {u) 



(E4) 



u + {nir/L) 2 KuM(u) 

Then, we have the solution in the following form: 
P(x,t) - T,n=o X n{x)T n (t). Finally, T„(0) in Eq. (E4) 
can be determined from the initial condition [Eq. (E2)] 
as 

T (0) = -, T„(0) = lC os-^ (n = l,2,...). (E5) 
Finally, we obtain 

, s 1 2 ,22, cos ^ cos ^ 

P(x,u) = — + -Y ^ 4 . E6 

Lu L £^ u + (nu: / L) 2 KuM(u) 

Note that P(x, u) in the above equation is just the tran- 
sition probability P(x 7 u;x s ,0) used in Sec. III. 



Appendix F: Riemann zeta function 

First, we define the functions C e (/?) an d as 

c e (/5)-E^> c°(/3)-E^- ( F1 ) 



n=2 
mcvcn 



n=l 
n:odd 



By using the integral representation of the gamma func- 
tion, 



1 1 



m Jo 



, dvv^~ 1 e~ nv 
" r(/3) ' 

it is easy to obtain the following formula: 

COS) 



2 /3 ' 



(F2) 



(F3) 



where is the Riemann zeta function. From Eq. (F3), 
we have 

2^-1 



C°(/3) 



2/3 



-C(/?)- 



(F4) 



In particular, we have C°(4) = tt 4 /96 and C°(2) = tt 2 /8. 



Appendix G: Ergodic theorems for general 
observables 



In this appendix, we summarize the ergodic properties 
of general observables. As an example, we derive the 
spatial distribution of free CTRWs. First, we consider a 
renewal process with renewal times tk (ft = 1, 2, ... ) (see 
Sec. II). We also define observables that take nonzero 
values only at the renewal times t' = tk (k = 1, 2, . . . ): 



h(t') = j2 H kHt'-t k ), 



(Gl) 



k=l 
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where {Hk} are random variables with the same mean 
value (H k ) = [i h (k = 1,2, ...). {H k } are assumed to be 
independent of N t , but not necessarily mutually indepen- 
dent. We also assume the ergodicity with respect to the 
operational time, i.e., the number of renewals n: 



1 ™ 

~y^ H k-^h, as 



oo. 



(G2) 



k=i 



This relation is just the law of large numbers. A suffi- 
cient condition for Eq. (G2) is that the correlation func- 
tion C(n) = {HkHk+n) — (Hk) (Hk+n) decays faster than 
n-T (7 > 0) [11, 65]. 

From Eq. (Gl), the time average of the function h(t') 
is given by 



1 



t i N t 

. , dt'h{t') = -Y j H k 

1 J 1 fc=i 



(G3) 



Note that the value of the RHS of Eq. (G3) depends on 
the trajectories of the renewal process in general. For 
example, if H k = 1 (k = 1,2,...) and the PDF of the 
renewal time r is give by a power law w(t) <~ 1/ r 1+a with 
< a < 1, the RHS of the above equation is equivalent 
to Nt/t, which follows the ML distribution [Eq. (41) with 
A = 0] as t — > oo. 

First, we consider the case in which /z/, ^ 0. In this 
case, we obtain the following equation from Eq. (G2): 



1 Nt 
— TH k 



— > (jbh, as t — > oo. 



fc=i 



Thus, Eq. (G3) can be rewritten as 

dt h(t ) ~ —Hh 



t 



t 



(G4) 



(G5) 



Similarly, if another observable g(t') defined by g(t') = 
J2T=i Gk&{t' — t k ) satisfies the same conditions as h(t'), 
we have 

I^W^f^ as ^ (G6) 

/„* h(t>)dt> Hh 

Note that the RHS is not a random variable. Equation 
(G6) is a stochastic version of Hopf's ergodic theorem for 
dynamical systems [28]. 

Next, we consider the case \ih = 0. We assume that the 
correlation function C(n) decays faster than 1/n; then, 
from the central limit theorem, we have 



1 " 



as n — > oo, 



(G7) 



k=i 



where F a is a random variable following the Gaussian 
distribution with mean and variance a 2 [11, 66]. The 



variance a 2 is given by a 2 = C(0) + 2^ =1 C(n) [11]. 
From Eqs. (G3) and (G7), we obtain 



dt'h(t') 



1/2 



-Fa- 



(G8) 



t Jo t 
If we define a random variable Y t as Y t = (N t /t a ) 1 F a , 

its characteristic function (e^ y '} is given by 



a 2 e 



(G9) 



where f\(£,t) is the Laplace transform of the PDF 
f\(x,t) [Eq. (41)]. Although, in general, it is difficult 
to obtain the PDF of Y t , we can derive it explicitly for 
A = 0. In this case, since fo(0 = /o(£>i)i which is the 
Laplace transform of the ML distribution, we obtain 



/o(o = E 



k=0 



1 



c) r(afe + l)' 



(G10) 



To derive this equation, let us introduce an auxil- 
iary variable h in Eq. (38) as Prob (N t /t a < xh) = 
ffa-i/a dr Ptl ( t ,0). Then double Laplace transforma- 
tions with respect to h and x give X^feLo ( — £/ csQ ) fe j 
where s is the Laplace variable conjugate to h. Finally, 
taking the inverse Laplace transform with respect to h 
and setting h = 1, we obtain Eq. (G10). Here note that 
Prob (N t /t a < xh) is the integral of the PDF f (x, t). By 
using Eqs. (G9) and (G10), we obtain the characteristic 
function for A — explicitly: 



OO / 



2A- 



1 



f^W^c/aJ T((a/2).2k + lY 



(Gil) 



Note that the PDF of Y t is symmetric with respect to 
the y-axis from its definition and that the even order 
moments of Y t are equivalent to those of the ML distribu- 
tion with index a/2 and scale factor \[2cja [Eq. (G10)]. 
This means that the PDF of Y t , fo(x), is a symmetric 
extension of the ML distribution; therefore we obtain 



1 00 F ( ka 
k=l 



1 



2c 



ife-i 



fc! 



sin 



(G12) 

for — oo < x < oo. The PDF /q(x) is called the bilateral 
Mittag-Leffler distribution [52]. The simplest example of 
this PDF is the spatial distribution of one-dimensional 
free CTRWs. Let the random variable H k be the dis- 
placement of the fc-th jump (e.g., H k = ±1 for CTRWs 
with jumps only to the nearest neighbor sites). If the 
jumps are symmetric, then (Hk) = 0. Thus, the time 
average [Eq. (G8)] becomes a rescaled spatial position 
and the PDF [Eq. (G12)] is the spatial distribution of 
the free CTRWs. Note that this PDF is equivalent to 
the one derived from the FFPE [Eq. (46) in [64]]. 
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